Solutions to the Jaynes-Cummings model without the rotating-wave approximation 
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By using extended bosonic coherent states, the solution to the Jaynes-Cummings model without 
the rotating-wave approximation can be mapped to that of a polynomial equation with a single 
variable. The solutions to this polynomial equation can give all eigenvalues and eigenfunctions of 
this model with all values of the coupling strength and the detuning exactly, which can be readily 
applied to recent circuit quantum electrodynamic systems operating in the ultra-strong coupling 
regime. 
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The Jaynes-Cummings (JC) model [l| describes the in- 
teraction of a two-level atom (qubit) with a single bosonic 
mode. It is a fundamental one in quantum optics. Based 
on the assumption of nearresonance and relatively weak 
atom-cavity coupling, the rotating-wave approximation 
(RWA) is usually employed, and analytically exact solu- 
tion can be trivially obtained. 

Recently, the JC model is closely related to condensed 
matter physics. It can be realized in some solid-state 
systems recently, such as one Josephson charge qubit 
coupling to an electromagnetic resonator [2j, the super- 
conducting quantum interference device coupled with a 
nanomechanical resonator 0, 0], and the most recently 
LC resonator magnetically coupled to a superconduct- 
ing qubit 0-0] ■ In traditional quantum optics where the 
coupling between the two- level "natural" atom and the 
single bosonic mode is quite weak, RWA is the most 
useful approximation. However, in the circuit quantum 
electrodynamic (QED), the artificial atoms may interact 
very strongly with on-chip resonant circuits |6l4l0| . the 
RWA can not describe well the strong coupling regime @. 
Therefore, the JC model without the RWA is the focus 
of current interests Illl-l2ll . 



However, due to the inclusion of the counter-rotating 
terms, the Bosonic number is not conserved, the Bosonic 
Fock space has infinite dimensions, so any solution with- 
out the RWA is highly nontrivial. In the recent years, 
several non-RWA approaches have been proposed in the 
Dicke model[l2j], the quantum Zeno effect (2^]. and the 
spin-boson model[23j]. Especially, by using extended 
bosonic coherent states, the present authors have solved 
the Dicke model without the RWA exactly in the numeri- 
cal sensepjl- The most simple N = 1 Dicke model is just 
the JC model. This numerically exact solutions to the JC 
model are also described in detail in Refs. [H, [l9j]. To 
the best of our knowledge, an analytical exact solutions 
are still lacking in the literature to date. 

In this paper, we propose a new method to solve ex- 
actly the JC model without the RWA by means of ex- 



tended bosonic coherent states. The correlations among 
bosons are added step by step until further corrections 
will not change the results. Different from our previ- 
ous work where the pure coherent state are fixed, the 
eigenvalue a of the pure coherent state in this paper is 
tunable. By solving Schrbdinger equation, the expanded 
coefficients can be expressed by a through a recurrence 
relation. We then derived a polynomial equation with 
only single variable a. The solutions to this polynomial 
equation can give exactly all eigenfunctions and eigen- 
values of the JC model without the RWA with arbitrary 
parameters. 

Without the RWA, the Hamiltonian of a qubit inter- 
acting with a single bosonic mode reads (ft = 1) 



(1) 



where a + and a are the bosonic annihilation and creation 
operators of the cavity, A and u are the frequencies of 
the atom and cavity, g is the atom-cavity coupling con- 
stant, and (Jk(k = x, y, z) is the Pauli matrix of the two- 
level atoms. For convenience, we can write a transformed 
Hamiltonian with a rotation around an y axis by an angle 



H = -^a x +a)a + g{a) +a)a z , (2) 

in the unit of u) = 1. Associated with this Hamiltonian 
is a conserved parity n, such that [H, H] = 0, which is 
given by 



n 



a x exp 



(3) 



where N = a + a is the bosonic number operator, n has 
two eigenvalues ±1, depending on whether the excitation 
number is even or odd. So the system has the correspond- 
ing even or odd parity. 
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We propose the following ansatz for the wavefunction 
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. (4) 



where c n is the expansion coefficient and a is the eigen- 
value of the coherent state and will be determined later. 
ty+ ( v &-) is the eigenfunction of the even(odd) parity 
with the eigenvalue +1(— 1). The Schrbdinger equation 
then gives 



M 



[a + a + g (a + + a)] c n (a + )" exp (aa + ) |0) 

ra=0 

exp (— aa + ) |0) 



A M 

n=0 



M 

= E ± J2cn(a + ) n exp(aa+)\0). (5) 

n=Q 

By using [a, a + ] = 1, equating the coefficients of the 
terms of (a + ) m exp (aa + ) |0) on both sides, the above 
equations for both even and odd parity can independently 
give the following identities for any to 



(to + ag) c m + (a + g) c m _i + (to + 1) c m +i 



(2a) 



j=o J 



E ± c„ 



(6) 



By careful inspection of Eq. (j4]) , one can find that cq is 
flexible in the Schrbdinger equation where the normaliza- 
tion for the eigenfunction is not necessary, so we select 
Co = 1.0. The linear term in a + in the Fock space can 
be also determined by the value of a in the pure coher- 
ent state exp(aa+) |0). It is useless that both ci and a 
contribute to the linear term, so we can set c\ — 0. Once 
the first two terms are fixed, the coefficients of the other 
terms higher than a + should be determined by solving 
Schrodinger equation. The constant term yields 



+ A 
E ± = ag ^—. 



(7) 



Inserting Eq. (JT)) into Eq. ([6]), we have the following 
recurrence equation 
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For to = M , the terms higher then (a + ) M exp (aa 
are neglected, we may set cm+i = 0, then we have 

M± y) Cm + c «- 1 



(9) 



Note that all coefficients c's can be expressed by one vari- 
able a through Eq. (jHJ, so this is a one- variable polyno- 
mial equation of degree M . This is the central equation 
in this paper, which roots would give the exact solutions 
to the JC model without the RWA. 

To obtain the true exact results, in principle, the trun- 
cated number M should be taken to infinity. Fortunately, 
it is not necessary. It is found that finite terms in state 
(l4| are sufficient to give exact results in the whole cou- 
pling range. Typically, the convergence is assumed to be 
achieved if the results are determined within very small 
relative errors when the truncated number M increases 
further. We like to stress here that increasing the value of 
M would almost not bring additional effort to solve Eq. 
(|9]) in an ordinary PC. The precision for the results is 
only limited to the machine accuracy for all cases tested. 

To have a sense about this method, we perform the 
first-order approximation (FOD) by considering M = 2. 
The only one coefficient to be determined is c%, which is 
easily obtained by the Eq. ([S]) ca = — ^ [(1 ± A) a + g]. 
The nonlinear equation then is given by 

T&ga 2 -(l±A)a-g = 

. The value of a for the even parity (+) and odd parity 
(— ) are 
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The corresponding energies are 
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^ = -^(l ± Jl-^UWi (12) 
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Note that the lower eigenvalue for the even parity (+), 
i.e. Eq. (|T2l ) with plus sign in the RHS, increases with 
g. It is physically unreasonable, because the ground- 
state(GS) energy should decrease with the qubit-cavity 
coupling. The corresponding coefficient C2 is considerably 
larger than Co = 1, indicating that the wavefunction is 
not converging. So this solution should be omitted. 

The GS energies in the FOD with different ratios of 
A/lu for M = 2, i.e. Eq. (fT2l ) with minus sign in the 
RHS, are plotted in Fig. [T]with solid circles. A general- 
ized RWA (GRWA) for JC model was performed recently 
[TTI ] and the derived expression for the energy levels is sig- 
nificantly more accurate than that with RWA. The GS 
energies by the GRWA with the same ratios of A/ui are 
also list Fig. Q] with solid curves for comparison. To 
show the accuracy, we also collect the GS energies by the 
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FIG. 1: (Color online) Comparisons of the GS energies by the 
present FOD (solid circles), GRWA (solid curves) , and ED 
(open circles) as a function of g at A/w = 0.5 (top), A/uj = 1 
(middle), and A/cj = 1.5 (bottom). 



numerically exact diagonalization (ED) in Bosonic Fock 

(a+) m 

states v . ; |0) with open circles. 



It is interesting that the GS energies by the present 
FOD are much more close to the numerical ED ones for 
the zero, positive and negative large detuning. The rel- 
ative difference for the GS energy is less than 10~ 3 for 
g < 0.5. One may deeply impressed by this good agree- 
ment, because the present FOD is only a preliminary ap- 
proximation. Note that the maximum coupling constant 
in the ultrastrong-coupling regime in the circuit QED re- 
ported recently jg, Q is around g = 0.12. When only 
the GS energy is concerned in these systems, the present 
very simple analytical expression Eq. (|12l ) with minus 
sign in the RHS should be very helpful. 

If we increase the truncated number M, we would ob- 
tain more accurate energy levels. On the basis of the 
Abel-Ruffini theorem that the general solution in radi- 
cals is impossible to polynomial equations of degree five 
or higher 24[ , the most accurate analytically expressions 
of the energy levels of the system might be obtained by 
set M = 4, which is however not shown here due to the 
complexity. Naively speaking, one could have 7 eigen- 
states for M = 4. Actually it is not that case. For the 
energy expression Eq. (J7J, only real roots for a in Eq. 
© is reasonable, some complex roots should be omitted. 

It is very crucial to obtain the real roots of Eq. © 
for sufficient large M where the general solutions do not 
exist. To achieve this goal, we plot a two dimensional 
diagram y = f(a), where f(a) is just the LHS of Eq. 
©. The schematic view of the solutions for the one- 
variable polynomial equation Eq. © for g — 0.1, 0.5, 
and 1 are presented in Fig. [3J The real roots are just the 
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FIG. 2: (Color online) Schematic view of the solutions for the 
one- variable polynomial equation Eq. @ for g — 0.1 (top 
panel), 0.5 (middle panel), and 1.0 (bottom panel). The left 
panel is for the even parity and the right for the odd parity. 
The truncation numbers are M — 59 for the even parity and 
M — 60 for odd parity. The insets show the enlarged view. 



crossing points of the curve y — f{a) and the straight 
line y = 0. 

It is observed universally that if M is a even (odd) 
number for the even (odd) parity, the lowest eigenvalue 
for the even (odd) parity is physically unreasonable, be- 
cause the expansion coefficients in the wavefunction (U) 
are not converging. This characteristic has shown up in 
the case of M = 2 discussed above. Except this lowest 
eigenvalue, it will be confirmed later that all other ob- 
tained eigenvalues of this model are true eigenvalues of 
model. So for a given model parameters g and A, we 
choose M to be a odd (even) number for the even (odd) 
parity so that the physically unreasonable solution does 
not appear. 

In Fig. [21 we choose M = 59 for the even parity and 60 
for the odd parity, which is sufficient large to ensure the 
convergence. The value of f(a) is very large in the small 
a regime, and decreases quickly as a increases. The real 
roots are clearly exhibited in this two-dimensional plot 
and can be obtained easily by Maple in the practical cal- 
culation. The number of the real roots are considerably 
less than M. E.g. for g = 0.1, we obtain totally 27 real 
roots with 13 roots for the even parity and 14 roots for the 
odd parity, which are corresponding to 27 energy levels. 
Note that the roots for a for the same energy levels gen- 
erally decrease with the coupling strength g. According 
to the energy levels Eq. ([7]), this is physically reasonable, 
because the total energy of this system should decrease 
with the coupling strength g. 
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FIG. 3: (Color online) The absolute value of the coefficients 
c n for g — 0.1 (bottom panel), g = 0.5 (middle panel), and 
g = 1 (top panel). The left panel is for the even parity and 
right for odd parity. The truncation numbers are M — 59 for 
the even parity and 60 for the odd parity. 



Once the roots are at hand, one natural question is 
whether the energy levels and the wavefunctions corre- 
sponding to these roots are really the exact ones to the 
JC model. Actually, for every model parameters, we in- 
crease M by 2 in each step until the values of the roots 
for a are not modified. Our criterion is that, if the rel- 
ative difference for a is less than 10~ 8 , we think M is 
large enough to give the exact solutions to the JC model. 

The exact solutions are also confirmed by the conver- 
gence of the expanded coefficients in the wavefunction 
(j4|), which are displayed in Fig. [3l where the absolute 
value of coefficients c n normalized to the maximum value 
of {c} for g = 0.1, 0.5, and 1.0 for the GS state and the 
first 5 excited states are plotted. It is very interesting 
that the coefficients c„ vanish after n — 40 for g = 0.1, 
20 for g = 0.5, and 10 for g = 1. It follows that M = 40 
for g = 0.1, M = 20 for g = 0.5, and M = 10 for g = l 
are large enough to get the exact solutions for these 6 
eigen states. Further increase of M could not change the 
wavefunction at all. In the practical calculations, if we 
choose the truncation numbers M = 59 for even parity 
and 60 for odd parity, we can get exact solutions for the 
first more than 20 eigen states for g > 0.1, which might be 
practically very useful when applied to the circuit QED 
system 0, Hi- 
lt is also interesting to note from Fig. [3] that the 
necessary truncation number M decrease as the cou- 
pling strength increases, indicating that the present ap- 
proach can be applied to the JC model in the arbitrary 
ultra-strong coupling regime. In sharp contrast with the 
present approach, in the numerically ED in the bosonic 
Fock states, the dimension the truncated subspace in- 



creases considerably with the coupling constant g, due to 
more photons are involved. In the present approach, a in- 
finite photons have been already included in the bosonic 
coherent states in the wavefunction (Q}. 

In the strong coupling limit, the first term in the JC 
Hamiltonian ^ can be neglected, so we have 



(14) 



A = a + g,B = a — g 



(15) 



we can observe that the linear term for the bosonic op- 
erator is removed, and only the number operators A + A 
and B + B are left. Therefore we can readily obtain the 
eigenfunctions 



with 



| to) 



(«+ + .gM r 



| TO) 



(16) 



10) 



(a+ - g/ujf' 



A > I 



— I"/B ' 

(17) 

and the eigenvalues = to — g 1 for the to state. Note 
that the eigenstates are twofold degenerate in the strong 
coupling limit. Comparing with the present wavefunction 
(j4]), we can find that the truncation number is just fixed 
to be to for the mth excited state. In the strong cou- 
pling regime, as the coupling strength g increases, only a 
fewer terms are needed to add to the TOth wavefunction 
Eq. (fl6|) for the TOth state. The stronger the coupling 
strength is, the fewer additional terms arc needed, which 
highlights the advantage of this method. 

Finally, to show the exact nature of the solutions di- 
rectly, we compare the present analytically exact solu- 
tions with those by numerical ED in Bosonic Fock states 
in Tables I, II and III. We only list the energies of the 
GS and the first 8 excited states for three typical detun- 
ings A/uj = 1,0.5, and 1.5 at different coupling strength 
where the RWA is usually invalid. It is very important to 
observe that the present results are the same as those by 
ED in all cases. We should pointed out here that in the 
numerical ED, the dimension of the truncated subspace 
is up to 10 5 at strong coupling, more computer mem- 
ory is needed to store the matrix elements, and several 
hours of CPU time is taken to perform the diagonaliza- 
tion. In the present approach, the dimension of the trun- 
cated subspace is only around 60, and the results are ob- 
tained within one second in any PC. More importantly, 
we only need to solve one polynomial equation with a 
single- variable. To the best of our knowledge, this is the 
first analytically exact solution to the JC model without 
the RWA for arbitrary coupling strengths and detunings, 
which is very simple and can be easily employed. 
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g=o.i 


g=0.5 


g=0.5 


g=i.o 


g=1.0 


E 


-0.754009629 


-0.754009629 


-0.856475589 


-0.856475589 


-1.27755156 


-1.27755156 


E 1 


0.223060923 


0.223060923 


-0.203235835 


-0.203235837 


-1.06857920 


-1.06857920 




0.768900004 


0.768900004 


0.600664430 


0.600664428 


-0.357181745 


-0.357181745 


Ez 


1.20197785 


1.20197785 


0.959669469 


0.959669469 


0.130321638 


0.130321638 


Ei 


1.78994281 


1.78994281 


1.47931569 


1.47931569 


0.871878630 


0.871878627 


E 5 


2.18237074 


2.18237074 


2.04745904 


2.04745904 


1.19903842 


1.19903842 


Ee, 


2.80950833 


2.80950833 


2.45896800 


2.45896801 


1.76995826 


1.76995826 


E 7 


3.16398011 


3.16398011 


3.00884499 


3.00884499 


2.18522105 


2.18522105 


Eg 


3.82785600 


3.82785600 


3.53494142 


3.53494142 


2.95127586 


2.95127586 



